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The  goals  of  this  project  were  to  (1)  analyze  numerical  phenomena  such  as  locking  and  boundary 
layers  occurring  in  the  modeling  of  elastic  bodies,  and  obtain  methods  with  robust  performance, 
(2)  extend  this  analysis  to  hierarchies  of  models,  and  (3)  continue  investigation  into  the  p  and 
h  —  p  FEM.  Specifically,  the  locking  of  a  hierarchy  of  plate  models  was  analyzed  to  show  that 
only  the  the  lowest  order  Reissner-Mindlin  model  effects  were  significant.  Essentially  locking-free 
h  —  p  mixed  methods  were  established  for  the  elasticity  problem,  Stokes  flow,  Reissner-Mindlin 
plate  model  and  Naghdi  shell.  The  h  —  p  FE  approximation  of  boundary  layers  was  analyzed. 
Optimal  convergence  estimates  for  the  3-d  p  version  boundary  element  method  were  obtained. 
Numerical  quadrature  in  the  p  version  was  analyzed  and  exponential  convergence  of  an  7i  —  p 
quadrature  scheme  for  singular  integrals  arising  in  boundary  element  and  vortex  methods  was 
established.  Wavelet  based  Galerkin  boundary  element  methods  as  well  as  a  convergent  FEM 
for  a  class  of  nonconvex  variational  problems  were  developed. 
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Introduction 

Problems  in  stress  analysis  often  involve  structures  that  can  be  decom¬ 
posed  into  canonical  thin  subdomains,  such  as  (possibly  laminated)  beams, 
plates,  shells,  etc.,  linked  by  junctions.  A  necessary  condition  for  good  per¬ 
formance  of  the  FEM  on  such  structures  is  its  robust  performance  on  each 
thin  component,  considered  separately. 

Our  primary  objective  has  been  to  fully  understand  how  to  model  the 
canonical  subdomains  and  how  to  design  FEM  discretizations  which  are  ro¬ 
bust  with  respect  to,  for  example,  the  thickness  parameter.  Such  discretiza¬ 
tions  must  avoid  problems  like  locking  and  the  existence  of  boundary  layers. 
Our  investigation  has  focussed  on  both  these  areas. 

A  second  objective  has  been  to  investigate  the  feasibility  of  adaptive, 
hierarchical  modeling  of  thin,  three-dimensional  domains.  By  this,  we  mean 
the  automatic  selection  of  a  model  from  a  hierarchy  of  two-dimensional  mod¬ 
els  that  converge  to  the  actual  solution  as  their  order  increases.  (This  is  in 
contrast  to  the  usual  egineering  technique  of  using  a  fixed  model,  such  as  the 
Reissner-Mindlin  model.)  Under  this  grant,  we  have  developed  computable 
a-posteriori  error  estimators  for  the  modeling  error.  Existing  p  version  FEM 
capabilities  then  provide  a  natural  implementation  of  the  hierarchical  mod¬ 
els.  This  work  has  applications  in  the  use  of  hierarchical  modeling  for  mul¬ 
tistructures,  where  localized  areas  (such  as  joints)  are  replaced  by  models  of 
higher  order  (or  even  the  3-d  problem),  while  models  of  minimal  complexity 
are  used  for  other  components. 

The  above  ties  in  well  with  our  continuing  work  on  various  aspects  of 
the  h  —  p  version.  In  this  regard,  we  have  also  investigated  the  p  version  of 
the  boundary  element  method  for  three-dimensional  problems,  mixed  h-p 
methods  for  the  Stokes,  elasticity,  Reissner-Mindlin  plate,  and  Naghdi  shell 
problems,  errors  in  the  p  version  due  to  numerical  quadrature,  and  h-p 
quadrature  for  singular  integrals. 

Moreover,  we  have  proposed  and  analyzed  a  wavelet  Galerkin  scheme  for 
the  discretization  of  boundary  integral  equations  of  the  first  kind  on  polygo¬ 
nal  domains  that  yields  a  “fast”  algorithm  with  superconvergence  properties 
for  the  solution  at  interior  points.  In  addition,  we  have  analyzed  and  im¬ 
plemented  a  finite  element  discretization  for  nonconvex  variational  problems. 

1.  Locking  and  Robustness  of  Standard  FEMs 

Locking  is  the  phenomenon  by  which  the  numerical  approximation  of 
parameter-dependent  problems  deteriorates  for  values  of  the  parameter  close 
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to  a  limiting  value. 

The  principal  investigators  M.  Suri  and  C.  Schwab,  in  joint  work  with  Ivo 
Babuska,  extended  previous  results  for  the  analysis  of  locking  phenomenon 
for  the  Reissner-Mindlin  plate  to  a  whole  hierarchy  of  plate  models.  It  was 
shown  in  [15]  that  there  is  no  additional  locking  present  for  plate  models  that 
are  of  higher  order  than  the  Reissner-Mindlin  plate.  Hence,  in  the  context 
of  hierarchical  modeling,  the  elements  only  have  to  be  designed  so  that  they 
properly  deal  with  the  locking  of  the  lowest  (Reissner-Mindlin)  case. 

Suppose  the  Reissner-Mindlin  plate  problem  is  approximated  over  a 
square  domain  with  a  uniform  triangular  or  rectangular  mesh,  using  the 
h  version  with  polynomials  of  degree  p  for  the  rotations  and  degree  q  for  the 
displacements.  Then  the  locking  and  robustness  orders  in  terms  of  the 
number  of  degrees  of  freedom,  is  summarized  in  Table  1,  from  [15]. 

Table  1:  Locking  and  robustness  for  the  h  version  with  uniform  meshes. 


Type  of 

Degree 

Degree 

Order  of  locking,  r 

Robustness  order,  1 

Element 

P 

q 

f{N)  =  0{Nn 

g{N)  =  0{N-^) 

Triangle 

1 

9  >  1 

r  =  1/2 

/  =  0 

(Vp) 

2  <  p  <  4 

q=p 

r  =  1 

1 

II 

P  >  5 

r  =  1/2 

;  =  (p  - 1)/2 

2  <p  <  3 

1  >p  +  l 

r  =  1/2 

/  =  (p  -  l)/2 

p  >  4 

1  >p  +  l 

r  =  0 

/  =p/2 

Product 

1 

1  >  1 

r  =  1/2 

1  =  0 

(Sp) 

p>2 

q>p 

r  =  l/2 

1  =  (P-  l)/2 

Trunk 

1 

9>  1 

r  =  1/2 

1  =  0 

(Q'p) 

2 

?  =  2,3 

r  =  1 

1  =  0 

q  >  4 

r  =  1/2 

1  =  1/2 

p  >  3 

q=p 

r  =  3/2 

l  =  {p-  3)/2 

q  >p  +  l 

r  =  1 

;  =  (p-2)/2 

In  a  second  work  on  locking  [18],  an  earlier  mathematical  definition  of 
locking  was  developed  into  precise  computable  ways  to  quantify  it.  The  dif¬ 
ference  between  h  and  p/h  —  p  techniques  in  combating  locking  was  brought 
out  using  various  computational  results  for  problems  of  engineering  interest 
such  as  nearly  incompressible  elasticity,  anisotropic  heat  flow,  and  problems 
over  thin  domains.  Various  issues,  such  as  the  difference  in  locking  when 
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different  quantities  of  interest  were  calculated,  the  increase  in  locking  when 
curved  elements  were  used,  and  the  superiority  of  high-order  methods  in 
combating  locking  were  discussed  in  a  computational  framework. 

2.  h  —  p  Approximation  of  Boundary  Layers 

Boundary  layers  are  an  important  phenomenon  in  the  solution  of  plate 
and  shell  problems.  Their  numerical  approximation  has  up  till  now  been 
treated  in  an  ad  hoc  fashion,  usually  by  (uniform)  over-refinement  near  the 
edges.  In  [19],  we  carried  out  a  comprehensive  investigation  on  how  h,  p  and 
h—p  refinements  could  be  best  designed  to  uniformly  capture  such  functions, 
for  any  thickness  d.  Our  work  is  particularly  useful  in  the  context  oi  h  —p 
version  codes,  such  as  MSC  PROBE  and  STRESSCHECK. 

Since  the  boundary  layer  behavior  manifests  itself  only  in  the  direction 
perpendicular  to  the  edge,  a  crucial  aspect  was  to  investigate  the  uniform 
approximation  of  1-d  boundary  layer  functions.  A  key  result  that  we  estab¬ 
lished  is  that  the  h  —  p  version  gives  uniform  exponential  convergence  if  the 
first  element  is  of  size  0{pd),  with  only  minimal  additional  elements  required 
to  mesh  the  remaining  domain.  Figures  1  and  2  illustrate  this  through  nu¬ 
merical  experiments  on  a  singularly  perturbed  second-order  elliptic  problem 
on  [-1,1]  with  exact  solution 

Ud{x)  =  1  -  cosh{x/d)/  cosh(l/d), 

which  has  boundary  layers  at  both  1  and  -1.  It  is  seen  that  for  different 
values  of  d,  the  3-element  mesh  outperforms  both  the  p  version  and  the 
optimal  h  version  mesh  (which  we  showed  must  be  exponentially  graded)  (it 
also  outperforms  the  h—p  version  with  more  elements). 

The  paper  [8]  was  presented  at  the  AMS  Mathematics  of  Computation 
conference  at  the  University  of  British  Columbia,  Canada,  August  9-13  by 
M.  Suri.  In  it,  we  showed  how  two-dimensional  elements  that  are  both  free 
of  locking  and  uniformly  approximate  the  boundary  layer  may  be  developed 
for  the  Reissner-Mindlin  plate.  This  work  has  been  extended  in  the  Ph.D. 
thesis  of  Christos  Xenophontos,  which  is  currently  under  preparation. 

The  paper  [20]  was  presented  at  the  ICOSAHOM  meeting  in  Houston, 
TX,  June  5-9.  In  it,  we  extended  the  work  in  [19]  to  a  model  one-dimensional 
advection-diffusion  problem,  showing  that  a  suitable  variational  formulation 
with  the  h—p  version  once  again  leads  to  exponential  convergence. 

The  technical  brief  [21]  discusses  mesh  design  for  structural  plates.  We 
showed  how  the  improper  mesh  design  in  programs  such  as  STRESSCHECK 
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can  lead  to  completely  incorrect  predicted  shear  forces,  due  to  unsatisfactory 
resolution  of  the  boundary  layer.  Following  our  theory  in  [19],  the  mesh 
should  be  refined  near  the  edges,  a  strategy  that  has  been  explicitly  adopted 
in  the  STRESSCHECK  manual. 

3.  Mixed  and  Reduced-Constraint  h  —  p  Methods 

Although  the  p  version  is  generally  free  of  locking,  various  commercial 
h-p  codes  now  offer  the  user  a  choice  of  both  h  and  p  refinement,  the  idea 
being  that  mesh  refinement  with  low  p  is  combined  with  low  refinement  and 
high  p,  depending  on  the  nature  of  the  subdomain.  It  is  therefore  essential 
that  the  method  used  be  free  of  locking  in  areas  of  low  p  as  well,  i.e,  where 
the  convergence  is  achieved  primarily  by  mesh-refinement  {h  version).  This 
suggests  a  modification  of  the  underlying  variational  method,  to  reduce  the 
effect  of  the  locking  constraint.  Such  a  modified  variational  form  is  often 
more  accurate  in  recovering  quantities  of  interest  such  as  the  pressure  and 
the  stress  as  well. 
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(a)  Elasticity  and  Stokes  Problem 

In  joint  work  with  R.  Stenberg  [12],  we  investigated  the  mixed  h  —  p  version 
for  the  problem  of  nearly  incompressible  elasticity,  and  its  limit,  the  Stokes 
flow  problem.  The  mixed  h  —  p  version  combines  the  power  of  both  high- 
order  and  mixed  methods  to  combat  locking.  Unlike  the  standard  method, 
where  only  locking-free  approximations  in  the  displacements  (or  velocities) 
can  be  guaranteed,  the  mixed  h  —  p  method  can  be  designed  to  control  the 
locking  in  the  stresses  (or  pressures)  as  well. 

We  analyzed  several  families  of  rectangular  mixed  elements  (including 
those  in  codes  such  as  “NEKTON”  and  “P3CFD”),  which  had  previously 
been  analyzed  either  in  terms  of  only  the  h  version  or  only  the  p  version,  but 
not  both.  We  showed  that  for  the  popular  spectral  element,  “Qn-Qn-2” 
the  /i-approximability  for  the  pressures  is  not  optimally  balanced  against 
that  for  the  velocities.  We  formulated  a  new  mixed  h  —  p  element  in  which 
this  balance  is  optimal,  and  which  is,  moreover,  minimal  in  the  number  of 
degrees  of  freedom  used,  while  still  being  locking-free  for  both  velocities  and 
pressures.  Our  elements  have  been  implemented  in  a  fully  adaptive  h—p 
method  for  fluid  flow  by  J.T.  Oden’s  group  at  TICAM,  with  excellent  results. 
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(b)  Reissner-Mindlin  Plate 

In  [14],  we  extended  our  results  above  to  obtain  locking- free  reduced  con¬ 
straint  h  —  p  elements  for  the  Reissner-Mindlin  plate.  Our  elements  are 
variants  of  the  well-known  MITC  elements,  which  had  previously  only  been 
shown  to  be  h  stable. 

(c)  Naghdi  Shell 

In  [13],  we  formulated  a  reduced-constraint  h  —  p  finite  element  method  on 
rectangular  elements  for  the  Naghdi  shell  model.  This  model  suflFers  from 
locking  in  the  bending-dominated  case.  We  showed  that  by  altering  the  un¬ 
derlying  variational  form,  a  method  that  was  substantially  free  of  locking 
both  in  terms  of  h  and  p  was  obtained. 

(d)  Design  Principles  for  Reducing  Plate  Constraints 

Joint  work  with  J.  Pitkaranta  [22]  was  motivated  by  the  locking  constraints 
found  in  shell  problems.  So  far,  satisfactory  mixed  method  techniques  to 
handle  such  constraints  have  not  been  developed.  A  natural  strategy  would 
be  to  extend  the  derivation  of  plate  elements  to  the  shell  setting.  However, 
mixed  methods  for  the  Reissner-Mindlin  plate  are  derived  by  appealing  to  a 
Hemholtz  decomposition,  an  indirect  technique  without  an  apparent  analog 
for  shells. 

In  our  work,  we  developed  a  new  means  of  deriving  mixed  elements  for 
the  Reissner-Mindlin  plate,  which  does  extend  to  shells.  For  the  first  time, 
the  reduction  operators  and  polynomial  spaces  used  to  treat  the  constraints 
are  directly  derived.  Our  analysis  gives  minimal  conditions  on  plate  elements 
needed  to  prevent  locking,  and  also  to  approximate  the  boundary  layers 
present.  Using  our  analysis,  a  clear  picture  emerges  for  the  comparison  of 
currently  available  low-order  {h  version)  plate  elements,  such  as  the  MITC 
and  Arnold-Falk  ones. 

We  have  now  started  applying  our  method  to  the  problem  of  shell  locking. 
4.  Hierarchic  Modeling 

In  the  area  of  hierarchic  modeling  of  thin  structures  our  principal  objec¬ 
tive  has  been  the  development  of  computable^  a-posteriori  estimators  of  the 
modeling  error,  i.e.  the  error  between  the  solution  of  the  three-dimensional 
problem  and  that  of  any  model  in  the  hierarchy. 
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(a)  Heat  Conduction  in  a  Thin  Plate 

C.  Schwab  has,  in  joint  work  with  1.  Babuska,  completed  an  analysis  of  hier¬ 
archical  models  for  the  heat  conduction  problem  in  a  thin  plate  [5].  There,  a- 
posteriori  estimators  for  the  modeling  error  were  obtained  and  their  asymp¬ 
totic  (i.e.  as  the  thickness  of  the  structure  tends  to  zero)  and  spectral  (i.e. 
as  the  spectral  order  of  the  model  tends  to  infinity)  exactness  were  proved. 
Further,  computable  a-posteriori  bounds  on  the  effectivity  indices  for  these 
estimators  were  obtained  for  the  first  time.  Using  symbolic  manipulation, 
all  bounds  in  our  theory  were  shown  to  be  sharp.  These  bounds  allow  the 
assessment  of  the  reliability  of  the  modeling  error  estimators  for  each  given 
boundary  value  problem  based  on  the  computed  solution.  Moreover,  these 
estimators,  based  on  the  residual  tractions  on  the  faces  of  the  plate,  were 
shown  to  be  locally  asymptotically  exacts  i.e.  their  local  size  gives  a  reliable 
indication  of  the  subdomains  in  which  the  plate  model  needs  improvement. 
This  allows  one  to  adaptively  select  different  model  orders  in  different  parts 
of  the  domain. 

The  computational  feasibility  and  efficiency  of  this  strategy  has  been 
demonstrated  in  the  paper  [6],  where  it  was  shown  that  the  estimators  in¬ 
dicate  reliably  the  boundary  and  interior  layers  present  in  the  three  dimen¬ 
sional  problem,  and  that  the  local  increase  of  the  model  order  is  computa¬ 
tionally  feasible  and  leads  to  near  optimal  reduction  of  the  modeling  error. 
It  was  found  essential  for  the  reliability  of  these  modeling  error  estimators 
that  an  accurate  numerical  solution  of  the  currently  adopted  two  dimensional 
model  be  available,  underlining  once  again  the  necessity  of  investigating  the 
locking  and  boundary  layer  phenomena  for  the  whole  hierarchy  of  the  mod¬ 
els  in  a  unified  fashion. 

(b)  Isotropic  Elasticity  and  Monoclinic  Plates 

The  above  results  were  extended  to  hierarchic  plate  models  of  uniform  model 
order  in  the  context  of  homogeneous,  isotropic  elasticity  in  the  paper  [4]  and 
later  to  general,  monoclinic  plates  in  [24].  Here,  once  again,  estimators  that 
are  guaranteed  upper  estimators  can  be  obtained.  The  constants  character¬ 
izing  the  estimators  depend  on  the  elastic  moduli.  We  provided  two  ways 
for  computing  them: 

a)  we  gave  analytic  expressions  for  the  constants  by  direct  asymptotic 
analysis  in  [24], 

b)  we  showed  how  to  compute  them  numerically  by  solving  (once  and 
for  all)  one  generalized  eigenvalue  problem  on  the  unit  cube. 

In  the  paper  [23]  we  also  proved  the  asymptotic  and  spectral  exactness 
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of  the  modeling  error  estimators  with  the  constants  obtained  according  to 
method  a).  We  expect  the  approach  b)  to  be  more  general  than  a)  and  to 
apply  also  for  variable  thickness  plates  and  shells. 

Below  is  an  important  result  for  homogeneous  and  isotropic  plates  Q.  — 
u)  X  (— t/2,t/2)  of  thickness  t  with  polynomial  body  force.  It  shows  how 
the  membrane  and  bending  modeling  error  in  the  energy  norm  may  be  esti¬ 
mated  by  a  computable  estimator.  Here,  n  =  (ni, n2,  n^)  is  the  model  order, 
denoting  the  polynomial  degrees  ni,  712,713  chosen  in  the  three  displacement 
directions  respectively. 

~  ^  \\e(Q)  ~  ‘^mem\\E{Q.)  W^hen  “  '^ben^E{Q.) 


Here  is  an  elastic  modulus  and  the  constants  A^en?  ^mem  depend  only 
on  the  model  order  n  and  the  Poisson  Ratio.  They  must  be  computed  by 
solving  (once  and  for  all)  a  generalized  eigenvalue  problem  on  the  unit  cube 
in  IR^.  The  are  the  computable  residual  tractions  on  the  faces  of  the 
plate. 

To  illustrate  the  performance  of  (1),  we  modeled  (using  PEGASYS)  the 
bending  of  a  clamped  steel  plate  with  a  central  hole  due  to  normal  loading 
(isotropic  material,  E  =  10^,  u  =  0.3,  thickness  t  =  0.1a,  a=edge  length). 
To  compute  the  solution  of  the’  n  models,  a  very  simple  mesh  in  was  used 
with  thin  elements  near  the  edge  included  for  boundary  layer  resolution. 
The  modeling  error  estimate  (1)  was  compared  to  the  true  one  (measured 
against  a  3  -  d,  h  -  p-FEM  solution).  Table  1  shows  estimated  and  true 
modeling  errors  and  the  effectivity  indices  ©  =Estimator/||'U*^  —  for 

models  up  to  order  6.  We  see  that  (1)  always  overestimates  the  modeling 
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Table  2:  Performance  of  the  residual  modeling  error  estimators  for  various 
model  orders. 


error  with  quite  moderate  efFectivity  indices  and  that  (1)  is  accurate  even 
for  higher  model  orders.  The  theory  has  been  generalized  in  the  meantime 
to  hierarchic  plate  models  with  a)  orthotropic  and  b)  laminated  materials, 
which  are  available  in  PEGASYS.  These  will  be  presented  in  a  monograph 
which  is  presently  under  completion. 

(c)  Fourier  Analysis  Techniques 

In  [2],  we  developed  the  tool  of  Fourier  analysis  of  the  modeling  error.  This 
turned  out  to  be  quite  useful  in  the  analysis  of  locking  phenomena  in  [15], 
and  has  been  extended  to  include  arbitrary  model  orders.  The  theory  has 
been  generalized  in  the  meantime  to  hierarchic  plate  models  with  general 
monoclinic  materials  ([24])  which  are  available  in  PEGASYS. 

(d)  Subspace  correction  methods 

In  [17],  we  started  work  on  the  convergence  analysis  of  subspace  correction 
methods  for  the  iterative  solution  of  hierarchic  plate  models.  There  the 
connection  between  the  selection  of  the  transverse  shape  functions  in  the 
hierarchic  plate  models  and  the  convergence  rate  of  certain  block  iterative 
procedures  was  analyzed,  again  with  the  tool  of  Fourier  analysis  developed 
in  [2].  It  was  shown  that  a  simple  orthogonalization  process  applied  to  the 
transverse  shape  functions  will  yield  iteration  procedures  with  a  convergence 
rate  which  improves  as  the  plate  gets  thinner.  In  recent  related  work  by  a 
group  in  the  U.K.,  it  is  actually  shown  that  this  selection  of  basis  functions 
also  ensures  a  convergence  rate  independent  of  the  model  order,  i.e.  we  have 
a  spectrally  uniform  convergence  rate.  For  the  near  future  (1995/6),  exten¬ 
sion  of  these  results  of  elasticity  problems  is  anticipated. 

(e)  Boundary  layers  in  hierarchic  models 

In  [16],  we  analyzed  in  detail  the  boundary  layer  structure  in  hierarchic 
models  of  homogeneous,  isotropic  plates.  This  was  done  to  obtain  insight 
into  the  discretization  of  the  plate  models  themselves  in  connection  with  our 
h  —  p  analysis  in  [19]. 

5.  The  p  Version  Boundary  Element  Method 

The  BEM  is  increasingly  being  used  for  approximating  partial  differen¬ 
tial  equations  over  polyhedral  domains.  Optimal  rates  of  convergence  for  the 
p  version  BEM  for  the  singularities  in  the  case  of  three-dimensional  prob¬ 
lems  were  established  in  the  paper  [25].  In  this  work,  we  described  the  three 
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types  of  singular  components  that  arise  in  the  solution:  edge,  vertex  and 
edge-vertex  singularities,  and  analyzed  the  approximation  of  the  traces  of 
each  of  these  by  polynomial  subspaces  on  the  boundary.  Asymptotic  rates 
of  convergence  as  p  — >  oo  that  were  optimal  in  the  norm  over  the  bound¬ 
ary  for  0  <  t  <  1  were  established  and  applied  to  a  model  Neumann  problem. 

6.  Quadrature  Results 

(a)  h  —  p  quadrature  of  singular  and  near  singular  integrals 

In  [7],  C.  Schwab  applied  ideas  from  h—p  finite  element  methods  to  numer¬ 
ical  evaluation  of  singular  and  near  singular  integrals  which  arise  in  many 
engineering  applications  (BEM,  vortex  methods,  etc.).  Exponential  conver¬ 
gence  with  respect  to  the  number  of  quadrature  points  was  shown  to  hold 
uniformly  for  integrands  in  the  countably  normed  space  containing 

in  particular  all  weakly  and  near  singular  integrals  arising  in  3  -  d  BEM. 

(b)  Error  in  p  version  in  presence  of  quadrature 

The  question  of  which  quadrature  rule  to  use  is  an  extremely  important  one 
in  practice,  since  in  commercial  programs,  the  major  portion  of  time  spent 
is  on  numerically  calculating  the  stiffness  matrices  (especially  in  three  di¬ 
mensions).  Various  topics  related  to  numerical  quadrature  for  the  p  version 
were  investigated  by  Chang  Kim  as  part  of  his  Ph.D.  thesis  (under  M.  Suri), 
completed  in  Spring,  1992.  In  particular,  the  effect  of  using  curved  elements, 
unsmooth  data,  and  the  quadrature  error  in  lower  order  norms  were  inves¬ 
tigated  theoretically  and  computationally.  The  results  were  published  in  [1]. 

7.  Wavelet  Galerkin  Discretization  of  First  Kind  Boundary  Inte¬ 
gral  Equations  on  Polygons 

C.  Schwab  has,  in  joint  work  [3]  with  T.  von  Petersdorff,  shown  mathe¬ 
matically  and  computationally  that  an  optimal  order  convergent  discretiza¬ 
tion  of  first  kind  boundary  integral  equations  on  arbitrary  polygons  in  the 
plane  is  possible  in  0{N  log  N)  memory  and  0{N (log  iV)^)  operations  where 
N  denotes  the  number  of  degrees  of  freedom  on  the  boundary.  It  has  been 
shown  that  a)  the  corresponding  stiffness  matrices  have  bounded  condition 
numbers  in  wavelet  bases  and  b)  that  instead  of  iV^  nonzero  entries  only 
0{N  log  N)  nonzero  entries  need  be  computed  and  stored  while  the  optimal 
asymptotic  convergence  rate  is  preserved,  even  in  the  most  negative  norm 
on  the  boundary.  The  results  cover  interior  and  exterior  boundary  value 
problems  for  the  Poisson,  Helmholtz,  plane  elasticity  and  Stokes  equations 
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alike.  Spline  and  biorthogonal  wavelets  similar  to  those  of  Chui  and  Wang 
are  used.  The  results  have  been  extended  to  closed  surfaces  in  in  [23]. 

8.  FEM  for  Non-convex  Variational  Problems 

The  FEM  discretization  of  an  abstract  minimization  problem,  minF(u), 
related  to  Newton’s  problem  of  minimal  resistance  of  a  body  moving  through 
a  fluid  was  analyzed  in  [10].  The  functional  F  was  neither  convex  nor  grow¬ 
ing  at  oo.  It  was  shown  that  with  polygonal  domains  and  linear  Courant  tri¬ 
angles,  the  flnite  element  approximations  over  appropriate  admissible  classes 
converged  to  a  minimizer.A  class  of  projected  Newton  methods  for  the  dis¬ 
crete  problems  yielded  locally  superlinear  convergence.  Numerical  experi¬ 
ments  were  performed  for  a  model  functional  F. 

9.  Expository  Work  on  p  and  h  —  p  Versions 

A  joint  paper  written  by  I.  Babuska  and  M.  Suri  [9]  on  the  foundations 
of  these  methods  was  published  in  SIAM  Review.  The  paper  contains  sim¬ 
plified  proofs  and  numerous  illustrations  that  bring  out  the  main  theoretical 
and  numerical  characteristics  of  these  methods,  and  serves  as  an  introduc¬ 
tion  for  both  mathematicians  and  engineers. 

10.  Benefits  to  USAF  and  Technology  Transfer 

A  number  of  our  results  have  been  or  will  be  implemented  in  commercial 
h  —  p  codes  used  (among  others)  in  the  aerospace  industry. 

(a)  h  —  p  approximation  of  boundary  layers. 

The  mesh-degree  combinations  suggested  by  us  to  resolve  boundary  layers 
were  used  to  show  how  the  h  —  p  program  STRESSCHECK  (developed  by 
Prof.  B.  Szabo  at  ESRD,  Inc,  St.  Louis,  Mo.)  can  extract  quantities  of  in¬ 
terest  such  as  stresses  with  great  accuracy  even  in  the  presence  of  boundary 
layers.  ESRD  has  issued  a  technical  brief  based  on  this  work  [21]. 

(b)  Modeling  error  estimators. 

Our  results  in  the  papers  [4]  -  [6]  were  discussed  and  used  by  ESRD  Software, 
Inc.  in  St.  Louis,  MO  in  their  development  of  the  FE  code  PEGASYS  which 
features  advanced  modeling  capabilities  for  (laminated)  plates  and  shells.  It 
was  demonstrated  in  [“SBIR  Phase  I  final  report:  Hierarchic  Modeling  of 
laminated  plates  and  shells”,  R.  Actis,  ESRD  Software  Inc.,  St.  Louis,  Mo] 
that  the  estimators  proposed  in  [4]  can  be  easily  implemented  and  that  the 
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layers  and  singularities  require  enhanced  resolution  by  a  refined  model.  The 
capability  of  variable  model  order  and  its  self-adaptive  selection  analyzed  in 
[11]  will  be  implemented  and  also  generalized  to  nonlinear  problems  in  this 
commercial  environment  during  the  Phase  II  of  the  SBIR  project. 

(c)  Mixed  h-  p  methods. 

The  h  —  p  elements  developed  by  us  for  Stokes  problem  were  the  ones  chosen 
for  a  fully  adaptive  implementation  by  A.  Patra  and  J.  T.  Oden  at  TICAM. 

(d)  Modeling  of  multistructures. 

We  are  working  closely  with  MacNeal-Schwendler  on  this  project.  The  idea 
is  to  be  able  to  model  substructures  independently  and  then  join  them  to¬ 
gether  using  an  interface  technique.  The  end-product  will  be  an  /t  —  p  com¬ 
mercial  code  that  has  substructuring  capability. 

(e)  Numerical  quadrature  for  singular  BEM  integrals. 

Subroutines  based  on  the  results  on  numerical  quadrature  have  been  incor¬ 
porated  into  an  industrial,  3-d  boundary  element  code  called  TOLOPT’ 
of  ABB  (Asea  Brown  Boveri)  Research  division,  Heidelberg,  Germany,  for 
electrostatics  simulations. 

(f)  USAF  Ph.D.  student. 

Major  Lawrence  Chilton  of  the  USAF  is  doing  his  Ph.  D.  with  us.  He 
will  serve  on  the  mathematics  faculty  at  Wright-Patterson  AFB,  Ohio  upon 
completion. 


Participating  Professionals 

In  addition  to  the  Pis,  the  following  personnel  were  involved  in  the  research: 

(a)  Graduate  (Ph.D.)  Students 
Chang  Kim 

Christos  Xenophontos 
Major  Lawrence  Chilton,  USAF 
Padmanabhan  Seshaiyar 

(b)  Visiting  Mathematicians  and  Consultants 
Rolf  Stenberg,  Helsinki  Institute  of  Technology 
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(c)  Other  Collaborating  Mathematicians 
Ivo  Babuska,  University  of  Maryland  College  Park 
Juhani  Pitkaranta,  Helsinki  University  of  Technology 
Tobias  v.  Petersdorff,  University  of  Maryland  College  Park 
Bernd  Kawohl,  Universitaet  Cologne,  Germany 

Advanced  Degrees  Awarded 

Chang  Kim,  Ph.D,  UMBC,  1992.  Thesis  title:  “The  effect  of  quadrature  er¬ 
ror  in  the  p  version  of  the  finite  element  method.”  (The  results  are  available 
in  reference  [1].) 


Seminars  and  Conference  Talks 

Jun  ’92  ICOSAHOM  ’92  (International  Conference  on  Spectral  and  Higher 
Order  Methods),  Montpellier,  France  (June  21-25,  1992),  “The  p  ver¬ 
sion  of  the  boundary  element  method  over  polyhedral  domains”  (M. 
Suri)  (organizer  of  minisymposium,  “The  p  version  of  the  finite  ele¬ 
ment  method”) 

Nov  ’92  Brunei  University,  Uxbridge,  U.K.,  “Locking  and  robustness  in 
the  finite  element  approximation  of  elasticity  problems”  (M.  Suri) 

Jan  ’93  Oxford  University,  Oxford,  U.K.,  “Locking  and  robustness  in  the 
finite  element  approximation  of  elasticity  problems”  (M.  Suri) 

Apr  ’93  MAFELAP  7  (Seventh  International  Conference  on  the  Mathe¬ 
matics  of  Finite  Elements  and  their  Applications),  Uxbridge,  U.K. 
(April  27-30,  1993),  “The  numerical  resolution  of  boundary  layers  and 
locking  effects  in  the  Reissner-Mindlin  plate  model”  (M.  Suri)  (Pis 
co-organized  minisymposium,  “The  finite  element  approximation  of 
plates  and  shells” ) 

May  ’93  Helsinki  University  of  Technology,  Espoo,  Finland,  “Locking  ef¬ 
fects  in  the  finite  element  approximation  of  plate  models”  (M.  Suri) 

Jun  ’93  University  of  Jyvaskyla,  Jyvaskyla,  Finland,  “The  p  and  h-p  ver¬ 
sions  of  the  finite  element  method”  (M.  Suri) 
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Aug  ’93  Symposium  on  Mathematics  of  Computation  1943-1993:  A  Half- 
Century  of  Computational  Mathematics,  University  of  British  Columbia, 
Vancouver,  Canada,  August  9-13,  1993,  ‘‘Locking  and  boundary  layer 
effects  in  plate  models”  (M.  Suri) 

Aug  ’93  Second  U.S.  National  Congress  on  Computational  Mechanics,  Wash¬ 
ington,  D.C.  (August  16-18,  1993),  “Locking  effects  in  the  finite  ele¬ 
ment  solution  of  plate  models”  (M.  Suri) 

Aug  ’93  “The  optimal  convergence  rate  of  the  p-version  BEM  on  polyhe- 
dra”,  lABEM  ’93  conference  at  Braunschweig,  ERG.  (C.  Schwab) 

Nov  ’93  “A-posteriori  error  estimation  for  adaptive,  hierarchic  plate  mod¬ 
els”  (C.  Schwab),  invited  talk  at  the  Biannual  conference  on  compu¬ 
tational  mechanics  at  the  University  of  Hannover  (E.  Stein,  Chair) 

Mar  ’94  Georgia  Institute  of  Technology,  Atlanta,  Georgia,  “The  mixed 
h  —  p  method  for  problems  in  elasticity  and  Stokes  flow”  (M,  Suri) 

Mar  ’94  Brown  University,  Providence,  Rhode  Island,  “The  mixed  h  —  j) 
method  for  problems  in  elasticity  and  Stokes  flow”  (M.  Suri) 

Jul,  ’94  1994  SIAM  Annual  Meeting,  San  Diego,  CA  (July  25-29,  1994), 
“Mixed  h  —  p  finite  elements  for  elasticity  problems”  (M.  Suri) 

Oct  ’94  Boundary  Element  Methods:  Applications  and  Error  Analysis, 
Oberwolfach,  Germany  (October  2-8,  1994),  “The  optimal  p  version 
approximation  of  singularities  on  polyhedra  in  the  BEM”  (M.  Suri) 

Oct  ’94  Boundary  Element  Methods:  Applications  and  Error  Analysis, 
Oberwolfach,  Germany  (October  2-8,  1994),  “Wavelet  Methods  for 
Boundary  Integral  Equations”  (C.  Schwab) 

Mar  ’95  Texas  A  &  M  University,  College  Station,  TX,  “The  h  —  p  finite 
element  modeling  of  plates  and  shells”  (M.  Suri) 

Mar  ’95  University  of  Texas,  Austin,  TX,  “The  h  —p  finite  element  mod¬ 
eling  of  plates  and  shells”  (M.  Suri) 

April  ’94  (Conference  on  Multiscale  Methods  in  Numerical  Analysis  or¬ 
ganized  by  the  Weierstrass  Institute,  Berlin)  “Fully  Discrete  Wavelet 
Galerkin  Boundary  Element  Methods”  (C.  Schwab) 
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Jul  ’95  ICIAM  ’95  (International  Congress  of  Industrial  and  Applied  Math¬ 
ematics),  “Mixed  h  —  p  methods  for  plates  and  shells”  (M.Suri)  (Pis 
co-organized  minisymposium  entitled  “Modeling  Thin  Structures  in 
Mechanics  —  Asymptotic  and  Finite  Element  Analysis”) 

Jul  ’95  ICIAM  ’95  (International  Congress  of  Industrial  and  Applied  Math¬ 
ematics),  “Fully  Discrete  Wavelet  Galerkin  Boundary  Element  Meth¬ 
ods”  (C.  Schwab) 
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